\documentclass[10pt,letterpaper]{article}
\usepackage[latin1]{inputenc}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{moreverb}
\usepackage{graphicx}
\usepackage{fancyheadings}

\graphicspath{{./Figuras/}}

%-----------------------------------------------------------------------
% Text dimensions 
%-----------------------------------------------------------------------
\setlength{\textheight}{22cm}
\addtolength{\textwidth}{56pt} % Aproxx. ~16cm
%
% Becuase I don't need marginal notes, then: 
%
\setlength{\marginparsep}{0pt}
\setlength{\marginparwidth}{0pt}
%
% Header is equal to \textwdith
% This is used when \marginparsep and \marginparwith are not zero
%
%\addtolength{\headwidth}{\marginparsep}
%\addtolength{\headwidth}{\marginparwidth}
%
% Header is equal to \textwdith
%
\setlength{\headwidth}{\textwidth}
%
% We have to substrac some points on even pages,
% because we add some points to \textwidth
%
\addtolength{\evensidemargin}{-56pt}
%\addtolength{\oddsidemargin}{-10pt}

\setlength{\voffset}{-1cm}


%-----------------------------------------------------------------------
% Style of page and headings
%-----------------------------------------------------------------------
\pagestyle{fancyplain}
\renewcommand{\sectionmark}[1]
	{\markright{\thesection\ #1}}
\lhead[\fancyplain{}{\bf\thepage}]
	{\fancyplain{}{\bfseries\rightmark}}
\rhead[\fancyplain{}{\bfseries\leftmark}]
	{\fancyplain{}{\bf\thepage}}
\cfoot{}
%\lfoot[\fancyplain{}{\bfseries\thepage}]{\fancyplain{}{}}
%\rfoot[\fancyplain{}{}]{\fancyplain{}{\bfseries\thepage}}


\author{Luis M. de la Cruz}
\title{TUNA by example}

\begin{document}

\maketitle

\section{Introduction}
The present document describes the programs of the \verb|Examples| directory of TUNA library.
The problems that are solved are explained briefly. The main parts of the different codes are explained 
in detail. There are seven directories, each one collect examples of the same kind. There are
several files in each directory: source code (\verb|.cpp|), input files, Makefile and in some cases
directories where the results of the execution are stored. A complete problem is solved in one 
single program, and this program uses the tools developed in TUNA. The idea of these examples
is to show how those tools can be used and combined to solve convection problems in 
box-shaped domains.

\strut

\noindent
\textbf{Compiling and running}

\noindent
Without lost of generality, suppose you are in the \verb|03Laplace| directory. If you
type:

\begin{verbatim}
% make
\end{verbatim}

\noindent
the files \verb|hoffman.cpp|, \verb|laplace1D.cpp| and \verb|laplace2D.cpp| will be compiled.
You can instead compile one single program typing, for example:

\begin{verbatim}
% make laplace1D
\end{verbatim}

\noindent
the executable file will have the same name as the source file, in this case \verb|laplace1D|.
To run this executable type:

\begin{verbatim}
% ./laplace1D
\end{verbatim}

\noindent
the program will read some parameters from the corresponding input file and will present
results on the screen. Some of the examples can be displayed with GNUPLOT or OpenDX, and
this is controlled with the compiling options: \verb|-DWITH_GNUPLOT| or \verb|-DWITH_DX|, take
a look to the Makefile to determine which option can be used.

\noindent
To clean all the compilation files, including the executable, just type:

\begin{verbatim}
% make clean
\end{verbatim}

Mostly of the constructions in TUNA require the precision and the dimension as parameters.
These data are given in the format: \verb|<T_prec, Dim>|, where \verb|T_prec| can be
float, double or long double, and \verb|Dim| can be 1, 2 or 3.

\section{Directory: 01StructMesh}
The file \verb|uniform_meshes.cpp|, is a program that shows how to use  the class \verb|StructuredMesh<>|
 which is intended to handle structured meshes in Cartesian coordinates. In this example the class is 
specialized with the \verb|Uniform<>| adaptor to create uniform meshes. 

The headers we need to include in order to be able to use the specialization \verb|StructuredMesh<>| 
and the adaptor \verb|Uniform<>|, are \verb|StructuredMesh.hpp| and \verb|Uniform.hpp| respectively.
Both are located in the directory \verb|Tuna/include/Meshes/|  and are implemented as C++ template 
classes.

The \verb|input| file  contains: 1 1 1 5 5 4 which are the lenghts in the axes $x$, $y$, and $z$ 
respectively, and the number of nodes in the same order. The meshes constructed with these data 
are sketched in figure \ref{fig:mallas}.

\begin{figure}[h!]
$$\includegraphics[width=10cm]{mallas.png}$$
\caption{Uniform structured meshes.}\label{fig:mallas}
\end{figure}

An extract of the output generated by this program is showed below:

\begin{footnotesize}
\begin{verbatim}
>----> Mesh inicialization 

 +----- TUNA info -----+
 | StructuredMesh<  Uniform < 1D > >
 +-----+
 | X-axis : length = 1; nodes = 5; volumes = 6; delta = 0.25
 +---------------------+
 
 +----- TUNA info -----+
 | StructuredMesh<  Uniform < 2D > >
 +-----+
 | X-axis : length = 1; nodes = 5; volumes = 6; delta = 0.25
 | Y-axis : length = 1; nodes = 5; volumes = 6; delta = 0.25
 +---------------------+

... MORE INFO ABOUT MESHES ...

\end{verbatim}
\end{footnotesize}

This output is generated by the next extract of code:

\begin{footnotesize}
\begin{listing}{1}
    std::cout << "\n\n >----> Mesh inicialization \n";
    StructuredMesh<Uniform<float, 1> > mesh_1D(length_x, num_nodes_x);
    mesh_1D.print();
    StructuredMesh<Uniform<float, 2> > mesh_2D(length_x, num_nodes_x,
					       length_y, num_nodes_y);
    mesh_2D.print();
\end{listing}
\end{footnotesize}

Lines 2 and 4 construct the objects \verb|mesh_1D| and \verb|mesh_2D| which
will contain all the information about the meshes (nodes, volumes, etc.). The function
\verb|print()| display all these information, lines 3 and 6.
The objects can be used later to solve a problem on the meshes generated.

\section{Directory: 02NonUniformMesh}

Program \verb|nonuniform_meshes.cpp| construct structured non uniform Cartesian 
meshes using the adaptor \verb|NonUniform<>|. Next extract of code shows how this
is done:
 
\begin{footnotesize}
\begin{listing}{1}
    StructuredMesh<NonUniform<double, 1> > mesh_1D(length_x, num_nodes_x);
    mesh_1D.calcPoints(X, MyFunction(lambda));
    mesh_1D.calcDeltas(X);
\end{listing}
\end{footnotesize}

In line 1 the object \verb|mesh_1D| is created to manage all the operations and
attributes of a non-uniform mesh. Lines 2 and 3 construct the points and the deltas
in $x$ direction, respectively. The functor \verb|MyFunction()| is used to implement
the point distribution function as follows:

\begin{footnotesize}
\begin{listing}{1}
   class MyFunction {      
       double lambda; // distribution factor 
     public:                                    
       MyFunction(double l) : lambda(l) { }   
  
       double operator() (double epsi) {      
          return ( (exp(lambda * epsi) - 1.0) / (exp(lambda) - 1.0) ); 
       }
   };
\end{listing}
\end{footnotesize}

This functor implements the function $f(\xi) =(e^{\lambda\xi} - 1)/(e^{\lambda} - 1)$ in
lines 6-8, where the \verb|operator()| is overloaded. If another function is
required, we just need to change line 7 and implement such function. The function is evaluated
in every $\xi$, which goes from 0 to 1, in increments of $1 / N -1 $, where $N$ is the number
of nodes on the particular axis direction.

The file \verb|input| contains the length of the domain in each axis, 8  4 2, the number of nodes,  
11  9 6, and the value for $\lambda$, in this case equal to 2.     
 
The results of this program are written in the \verb|Data| directory.  The figure \ref{fig:mallas1} 
were generated with OpenDX, using the \verb|mesh.net| script located in the same directory. 
If you use the option for GNUPLOT the graphs appear as the program is executing, you need 
to type $<$1$>$ and then $<$return$>$ each time a graph is presented in order to show the next.

\begin{figure}[h!]
\begin{center}$
\begin{array}{cc}
\includegraphics[width=7.0cm]{mesh1D2D.png} &
\includegraphics[width=6.0cm]{mesh3D.png} \\
\end{array}$
\end{center}
\caption{Non-uniform meshes in 1D, 2D and 3D.}\label{fig:mallas1}
\end{figure}

\section{Directory: 03Laplace}

The files in this directory solve the Laplace equation : $\nabla^2 \psi = 0$. 
In the codes $\psi \equiv T$.

\begin{description}
\item[laplace1D.cpp] : This code solves the Laplace equation in one dimension.
This equation is solved in the interval $ x \in [0,1] $.
Initially we have $ T = 10 \cos(3 \pi x) $ and the boundary conditions are: 
$ T=10 $ for $ x=0 $, and $ T=0 $ for $ x = 1 $. 
\item[laplace2D.cpp] : This code solves the Laplace equation in two dimensions.
This equation is solved in the unit square $ (x,y) \in [0,1]\times [0,4] $.
The boundary conditions are:
$ T=0 $ for $ x=0, 4$ and $ T=10 \sin (\frac{\pi}{2} x) $
for $ y=0, 1$.
\item[hoffman.cpp] : This code solves the Laplace equation in two dimensions.
 This equation is solved in the rectangle
 interval $ [0,1] \times [0,2] $. The boundary conditions are: $T = 100$ for
 $y = 2$, and $T= 0$ on the other walls. The exact solution of this problem is: 
\[ 
T(x,y) = 10\sin(\pi x) \frac{e^{\pi x}-e^{-\pi y}}{e^{\pi}-e^{-\pi}}
\] 
\end{description}

All these problems start defining a mesh, for example in \verb|laplace2D.cpp| this is done
as follows:

\begin{footnotesize}
\begin{verbatim}
    StructuredMesh<Uniform<float, 2> > mesh(length_x, num_nodes_x, length_y, num_nodes_y);
\end{verbatim}
\end{footnotesize}

The arrays to store the solution are declared as:

\begin{footnotesize}
\begin{verbatim}
    ScalarField2D T (mesh.getExtentVolumes());
    ScalarField2D Tn(mesh.getExtentNodes());
\end{verbatim}
\end{footnotesize}

The first declaration gets the number of volumes from the object \verb|mesh| and use this info
to create an array of compatible size with the mesh. Similarly, in the second declaration, the number of nodes are
obtained from the object \verb|mesh| and the \verb|Tn| array is created. This latter is used to store and
to plot the results. Some tools for interpolation from nodes to center of volumes and vice versa
are included in TUNA, for example:

\begin{footnotesize}
\begin{verbatim}
    NumUtils::interpolateToVolumes(T, Tn);
\end{verbatim}
\end{footnotesize}

\noindent
interpolate the data from \verb|Tn|, info on the nodes, to \verb|T|, info on the center of volumes.

Then, the arrays to store the linear system containing the coefficients of the discretization are created. 
We need a matrix and a vector:

\begin{footnotesize}
\begin{verbatim}
    DiagonalMatrix<float, 2> A(num_nodes_x, num_nodes_y); 
    ScalarField2D            b(num_nodes_x, num_nodes_y);
\end{verbatim}
\end{footnotesize}

The object \verb|A| stores the diagonals of the matrix, where appear terms different from zero. 
The finite volume method generate tridiagonal matrices in one dimension, five diagonals in
two dimensions and seven diagonals in three dimensions. The size of the matrix and the vector
is calculated automatically using information from the number of nodes in each direction.

With all these information we are able to create the equation to be solved, and this is done as
follows:

\begin{footnotesize}
\begin{listing}{1}
    ScalarEquation< LaplaceScheme<float, 2> > laplace(T, A, b, mesh.getDeltas());
    laplace.setDirichlet(TOP_WALL);
    laplace.setDirichlet(BOTTOM_WALL);
    laplace.setDirichlet(LEFT_WALL);
    laplace.setDirichlet(RIGHT_WALL);
    laplace.setGamma(1.0);
    laplace.print();
\end{listing}
\end{footnotesize}

In line 1, the object \verb|laplace| is created from the specialization \verb|ScalarEquation<>|, which in turn
is adapted using \verb|LaplaceScheme<>|. This adaptor contains the definitions to calculate the coefficients
of the discretization using a particular numerical scheme. To create the object, the next arguments are
required: the scalar field where the solution will be stored, the linear system (\verb|A| and \verb|b|) and the
mesh of the domain. Then some parameters of the equation are defined: the boundary conditions, lines 2-5, and
the diffusion coefficient, line 6. Finally, the function \verb|print()| is used to display information of the 
equation, line 7.

Once the equation is defined, it is possible to calculate the coefficients, to fill the linear system and to
solve it. This is done with:

\begin{footnotesize}
\begin{listing}{1}
    laplace.calcCoefficients();
    int iter = Solver::TDMA2DX(laplace, 0.001, 30);
    laplace.update();
    NumUtils::interpolateToNodes(Tn, T);
\end{listing}
\end{footnotesize}

In line 1 the coefficients are calculated, in line 2 the algorithm \verb|TDMA2DX| is used to solve the
linear system. The solution is stored internally in the object \verb|laplace|, so in order to get that
solution stored in the array \verb|T| we need to apply the function \verb|update()|, line 3. Finally, the
information stored in \verb|T|, which is on the center of the volumes, is interpolated to
the nodes using the function \verb|interpolateToNodes(Tn, T)|, line 4.
Results of the programs \verb|laplace2D.cpp| and \verb|hoffman.cpp| are shown in figure \ref{fig:laplace}

\begin{figure}[h!]
\begin{center}$
\begin{array}{cc}
\includegraphics[width=2.5in]{laplace2D_01.png} &
\includegraphics[width=2.5in]{laplace2D_02.png} \\ 
\mbox{(a)} & \mbox{(b)} \\
\includegraphics[width=2.5in]{HoffmanNum_01.png} &
\includegraphics[width=2.5in]{HoffmanErr_01.png} \\
\mbox{(c)} & \mbox{(d)} \\
\end{array}$
\end{center}
\caption{(a)-(b) Results of laplace2D.cpp. (c) Numerical solution of hoffman.cpp and (d)
numerical error with respect to the analytical solution.}\label{fig:laplace}
\end{figure}

\strut

\noindent
\textbf{Observation}

The advantage of using these tools from TUNA is that all the examples apply the same sequence of operations. 
All the tools can be used in 1, 2 and 3 dimensions. In the next examples we only 
highlight the particular behavior of each problem.

\section{Directory: 04HeatDiffusion}
This code solves the next unsteady equation
\[
\frac{\partial \psi}{\partial t} = \nabla^2 \psi .
\]
In the codes $\psi \equiv T$.

\begin{description}

\item[heat1.cpp] :
The equation is solved in the interval $ x \in [0,1] $. The initial condition is 
$ T = 10 \cos (3 \pi x)$, and the boundary conditions are $T=10$ for $x=0$, and 
$ T=-10 $ for $ x = 1 $.

\item[heat2.cpp] :
The equation is solved in the unit square $ (x,y) \in [0,1] \times [0,1]$. The initial condition
is $T = 0$. The boundary conditions are: $T = 100$ for $y = 2$, and $T= 0$ on the other 
walls.

\item[heat3.cpp] :
This code solves the equation in $(x,y,z) \in [0,1]\times [0,1] \times [0,1]$
The top wall is at $T = T_c$ while bottom wall is at $T = T_h$, and $T_c < T_h$.
Other walls are adiabatic ($\partial T / \partial n = 0$). The initial condition
is $T = 0$.

\item[oops.cpp] 
The equation is solved in a unit square $(x, y) \in [0,1] \times [0,1]$. 
The initial condition is "OOPS" (Object-Oriented Programming for Science) 
inside the square. The boundary conditions are $T = 0$ in all walls 
of the square.
\end{description}

To illustrate how these problems are solved with TUNA, we explain the code 
\verb|heat3.cpp|. The mesh, the scalar fields and the initial conditions are defined in 
three dimensions:

\begin{footnotesize}
\begin{verbatim}
    StructuredMesh<Uniform<double, 3> > mesh(length_x, num_nodes_x,
                                             length_y, num_nodes_y,
                                             length_z, num_nodes_z);
// ...
    ScalarField3D T(mesh.getExtentVolumes());
    ScalarField3D Tn(mesh.getExtentNodes());
//...    
    Range all = Range::all();
    T(all, T.lbound(secondDim), all) = 0.5;  // Bottom wall
    T(all, T.ubound(secondDim), all) = -0.5; // Top wall
    NumUtils::interpolateToNodes(Tn, T);
\end{verbatim}
\end{footnotesize}

Then we create the linear system, the equation to be solved along with some important
parameters and the boundary conditions:

\begin{footnotesize}
\begin{verbatim}
    DiagonalMatrix< double, 3 > A(num_nodes_x, num_nodes_y, num_nodes_z); 
    ScalarField3D               b(num_nodes_x, num_nodes_y, num_nodes_z);
//...
    ScalarEquation<CDS<double, 3 > > heat(T, A, b, mesh.getDeltas());
    heat.setDeltaTime(dt);
    heat.setGamma(1.0);
    heat.setDirichlet(TOP_WALL);
    heat.setDirichlet(BOTTOM_WALL);
    heat.setNeumann(LEFT_WALL);
    heat.setNeumann(RIGHT_WALL);
    heat.setNeumann(FRONT_WALL);
    heat.setNeumann(BACK_WALL);
    heat.print();
\end{verbatim}
\end{footnotesize}

Finally, a loop is initiated to solve the problem in time:

\begin{footnotesize}
\begin{verbatim}
    while ( (error > tolerance) && (iteration < max_iter) )
    {
        heat.calcCoefficients();
        Solver::TDMA3D(heat, tolerance, tdma_iter);
        error = NumUtils::calcErrorL2(T, heat.phi);
        residual = heat.calcResidual();
        heat.update();
    }
\end{verbatim}
\end{footnotesize}

The results of \verb|heat2.cpp|, \verb|heat3.cpp| and \verb|oops.cpp| are 
shown in figure \label{fig:heat}

\begin{figure}[h!]
\begin{center}$
\begin{array}{cc}
\includegraphics[width=2.5in]{heat2.png} &
\includegraphics[width=2.5in]{heat3.png} \\ 
\mbox{(a)} & \mbox{(b)} \\
\includegraphics[width=2.5in]{oops1.png} &
\includegraphics[width=2.5in]{oops2.png} \\
\mbox{(c)} & \mbox{(d)} \\
\end{array}$
\end{center}
\caption{(a) heat2.cpp, (b) heat3.cpp, (c) first step of oops.cpp and (d)
result of oops.cpp after 100 time steps.}\label{fig:heat}
\end{figure}


\section{Directory: 05ConvDiffForced}

The examples of this directory solve the equation:

\[ 
\frac{\partial \psi}{\partial t} + 
u_j \frac{\partial \psi}{\partial x} = \frac{\partial^2 \psi}{\partial x^2}
\]. 

In the codes $\psi \equiv T$.

The programs are:

\begin{description}

\item[convection1.cpp] :
The equation is solved in the domain $ x \in [0,2.5]$. 
The velocity is equal to 1. The boundary conditions are $T = 1$ for $x=0$
and $T = 0$ for $x=2.5$.

\item[convection2.cpp] :
The equation is solved in the unit square $ (x, y) \in [0,1] \times [0,1] $. 
The velocity is prescribed and has the form:
\[
u = -A \cos(\pi y) \sin(\pi \lambda x / L_x );
\]
\[
v = \frac{A \lambda}{L_x} \sin(\pi y) \cos(\pi \lambda x / L_x);
\]

The top wall is at $T = T_c$ while bottom wall is at $T = T_h$, and $T_c < T_h$.
Other walls are adiabatic ($\partial T / \partial n = 0$). The initial condition
is $T = 0$.

\item[convection3.cpp] :
This problem is similar to the previous one but in three dimensions. In this case
the third component of the velocity is $w = 0$.


\end{description}


\section{Directory: 06ConvDiff}

\begin{description}
\item[convdiff1.cpp] :
\item[convdiff2.cpp] :
\end{description}

\section{Directory: 07ConvDiffLES}

The example \verb|convdiffLES.cpp| 

\section{Final Remarks}

\end{document}